** Spatial approach, 3: Define Neighborhood ** 

		**********------------------------------**********
	* Contents: 
		* Build weight matrices for a given neighborhood definition
		* Fan corresponds to "V-Shaped" in final paper  
		* Input: cell-level data based on spatial_2_collapse_s 
		* Output: weight matrixes in .dta format 

capture set maxvar 30000, permanently // required to run file.  
	
capture log close 
clear all
set more off 

cd "R:\WSV2\TBu_AKe\Spatial_NEW"

capture mkdir Data
global store "R:\WSV2\TBu_AKe\Spatial_NEW\Data"

log using spatial_3_wmat_fan, replace 


cd "R:\WSV2\TBu_AKe\Spatial_NEW"


use "C:\Users\hy65byfe\Desktop\smerge_0712\spatial_2_collapse_s.dta" // data set with 4950 cells, created in spatial_2_*  


scalar p = 0 // if 0, straight distance without decay : x^0 = 1, so all are equal.
scalar sigma = 47*0.59 // sigma at comp. line 
scalar kappa = 51.7*0.59 // kappa at comp. line 
scalar D = 3 // width of increment 

***************************
**** reduce dimensions ****
***************************

keep if year == 2011 & ccode == 1 // reduce data in memory 
  
sort cell // now sorted: first bunching (1-50), then restricted (51-1015) then rest. 

gsort -bunching +cell 
gen bcell = _n if bunching == 1 // identify B

gsort -restricted +cell 
gen rcell = _n if restricted == 1 // identify R


**********************************************
*** define upper and lower bounds of cells ***
**********************************************
		
	** gen mid-point for E 
gen a_mid = a_bin 
gen e_mid = `=kappa' + `=sigma'*a_mid //  

	** gen upper and lower level of E at mid-point (e.g. 6kg, 6.5kg ...)
replace e_bin = k_bin + a_bin*`=sigma' // gravity center: upper edge of mid-point 
gen e_ml = k_bin + a_bin*`=sigma' - `=D' // lower edge of mid-point 


replace eei_bin = (e_bin/(51.7+47*a_bin))* 100 // assert that bins are consistent with eei boundaries. 

	** label as safety net 
label variable e_bin "E_i at midpoint"

label variable a_min "a_i left"
label variable a_max "a_i right"
label variable k_bin "kappa (std. E)"
label variable e_mid "E_m"


keep eei_bin a_* e_* k_bin s_bin eei_bin cell bcell rcell // cut 

sort cell // order is crucial 

***********************************
******* start loops ***************
***********************************

**********************************	 
*  3a: neighborhood slope / fan *
**********************************

***********************************
*** gradients from unit circle  ***
***********************************
	
	** sw 
scalar g0 = tan((225*(_pi/180)))
scalar g1 = tan((202.5*(_pi/180)))
	** horizontal 
scalar g2 = tan((180*(_pi/180)))
	** mixed (nw) 
scalar g3 = tan((157.5*(_pi/180)))
scalar g4 = tan((135*(_pi/180)))
scalar g5 = tan((112.5*(_pi/180)))
	** vertical 
scalar g6 = -tan((90*(_pi/180))) // has to be minus due to quadrant switch at 90. 'approach' from left on unit circle.
scalar g7 = tan((67.5*(_pi/180)))
scalar g8 = tan((45*(_pi/180)))

forvalues i = 0/8 {
	scalar delta`i'= (1/g`i')
}	

* scalar list 

**************************************
* 3b: define scalar for every cell ***
**************************************


sum cell, meanonly 
scalar cmax=`r(max)'
local c = `r(max)' // scalar for each cell.  
	
sum bcell, meanonly 
scalar bmax=`r(max)'

sum rcell, meanonly 
scalar rmax=`r(max)'

	*** scale to normalize slope of compliance line to 1. From formula: EEI = (E / k + ba ) * 100
scalar sd_a = 100/`=sigma'
scalar sd_e = 100


forvalues i = 1/`=cmax' {
	* 5-digit cell-id
	scalar cell`i' = 10000 + `i'  
	
	* capacity 
	quietly sum a_bin if cell == `i'
	local a = r(mean)
	scalar a`i' = `a'/sd_a
	
	quietly sum a_min if cell == `i'
	local a = r(mean)
	scalar a_min`i' = `a'/sd_a
	
	quietly sum a_max if cell == `i'
	local a = r(mean)
	scalar a_max`i' = `a'/sd_a
	
	* std.  
	quietly sum s_bin if cell == `i'
	local s = r(mean)
	scalar s`i' = `s'
	
	quietly sum k_bin if cell == `i'
	local s = r(mean)
	scalar k`i' = `s'
	
	* eei 
	quietly sum eei_bin if cell == `i'
	local x = r(mean)
	scalar x`i' = `x'
	
	* e [four courners and midpoint]
	quietly sum e_bin if cell == `i'
	local e = r(mean)
	scalar e`i' = `e'/sd_e
	
	quietly sum e_low_left if cell == `i'
	local e = r(mean)
	scalar e_ll`i' = `e'/sd_e
	
	quietly sum e_up_right if cell == `i'
	local e = r(mean)
	scalar e_ur`i' = `e'/sd_e
	
	quietly sum e_up_left if cell == `i'
	local e = r(mean)
	scalar e_ul`i' = `e'/sd_e
	
	quietly sum e_low_right if cell == `i'
	local e = r(mean)
	scalar e_lr`i' = `e'/sd_e
	
	quietly sum e_mid if cell == `i'
	local e = r(mean)
	scalar e_mid`i' = `e'/sd_e
}


***********************************
**** 3c: define neighborhood ******	 
***********************************
		 
****************************	 
	* inclusion conditions 
	
forvalues k = 1/`=bmax' {
     quietly sum cell if bcell==`k'
     local i = r(mean)
	  * i contains cellid of cell i in B 	 
          forvalues l = 1/`=cmax' {
          quietly sum cell if cell==`l'
          local j = r(mean)
		  * j contains cellid of cell j in R 
		  
		  scalar n_`i'_`j'  	= cond(a`j' == a`i' 									& x`j' > x`i' &  inrange(x`j', 59, 87),	1,0) // same size (rook)
		  scalar ne_`i'_`j' 	= cond(abs(e`i'-e`j') <= 6 	& a`j' <= a`i' 		   		& x`j' > x`i' &  inrange(x`j', 59, 87),	1,0)  // same e (+/-D)
		  
		  *** quadrants *** 
		  scalar nw_`i'_`j' 	= cond(a`j' <= a`i' 	 	& e`j' >= e`i'				& x`j' > x`i' &  inrange(x`j', 59, 87),	1,0) // nordwest (smaller j, H_o)
		  scalar no_`i'_`j' 	= cond(a`j' >= a`i'    		& e`j' >= e`i' 				& x`j' > x`i' &  inrange(x`j', 59, 87),	1,0) // nordost (bigger j)
		  scalar sw_`i'_`j' 	= cond(a`j' <= a`i' 	 	& e`j' <= e`i'				& x`j' > x`i' &  inrange(x`j', 59, 87),	1,0) // southwest (smaller and lower)
          }
     } 


******************************* 
***  4a: defining neighbors. ***
*******************************

forvalues k = 1/`=bmax' {
     quietly sum cell if cell==`k'
     local i = r(mean)
	  * i contains cellid of cell 1	 
          forvalues l = 1/`=cmax' {
          quietly sum cell if cell==`l'
          local j = r(mean)
		 * j contains cellid of cell 2
		 
			** distance between two cells in e-direction: negative if second cell consumes more ** 
		  scalar slope_`i'_`j' 	= ((e`j'-e`i')/(a`j'-a`i'))*(sd_a/sd_e) // TRAP: sign: +/- abd: missing when same size. 
			** convert to angle 
		  scalar angle_`i'_`j' 	= cond(slope_`i'_`j' != . , atan(slope_`i'_`j')*(180/_pi), 90) // TRAP: sign. transfrom slope to angle on unit circle.  
          }
     }
	 
	 
forvalues k = 1/`=bmax' {
     quietly sum cell if bcell==`k'
     local i = r(mean)
	  * i contains cellid of cell i in B 	 
          forvalues l = 1/`=cmax' {
          quietly sum cell if cell==`l'
          local j = r(mean)
		  * j contains cellid of cell j in R 
			 
			** 2 sections in SW quadrant 
		  scalar s1_`i'_`j' 	= cond(abs(angle_`i'_`j') >= 22.5 	& abs(angle_`i'_`j') <= 45		& sw_`i'_`j' == 1,	1,0)  // section 1
		  scalar s2_`i'_`j' 	= cond(abs(angle_`i'_`j') >= 0 		& abs(angle_`i'_`j') <= 22.5	& sw_`i'_`j' == 1,	1,0)  // section 2
		  
		    ** 4 sections in NW quadrant 
		  scalar s3_`i'_`j' 	= cond(abs(angle_`i'_`j') >= 0 		& abs(angle_`i'_`j') <= 22.5	& nw_`i'_`j' == 1,	1,0)  // section 3
		  scalar s4_`i'_`j' 	= cond(abs(angle_`i'_`j') >= 22.5 	& abs(angle_`i'_`j') <= 45		& nw_`i'_`j' == 1,	1,0)  // section 4
		  scalar s5_`i'_`j' 	= cond(abs(angle_`i'_`j') >= 45  	& abs(angle_`i'_`j') <= 67.5	& nw_`i'_`j' == 1,	1,0)  // section 5
		  scalar s6_`i'_`j' 	= cond(abs(angle_`i'_`j') >= 67.5	& abs(angle_`i'_`j') <= 90		& nw_`i'_`j' == 1,	1,0)  // section 6
		  
		    ** 2 sections in NO quadrant 
		  scalar s7_`i'_`j' 	= cond(abs(angle_`i'_`j') >= 67.5 	& abs(angle_`i'_`j') <= 90		& no_`i'_`j' == 1,	1,0)  // section 7
		  scalar s8_`i'_`j' 	= cond(abs(angle_`i'_`j') >= 45 	& abs(angle_`i'_`j') <= 67.5	& no_`i'_`j' == 1,	1,0)  // section 8
		 
		  
          }
     }	 	 
	 
	 

*********************************
*** 4b: transform to matrix *****
*********************************
display cmax // check dimension 

*** neighborhood matrices
*mat define X  = I(`=cmax') // knock out criterium
mat define S1 = I(`=cmax') // sw 1   
mat define S2 = I(`=cmax') // sw 2
mat define S3 = I(`=cmax') // nw 1 
mat define S4 = I(`=cmax') // nw 2
mat define S5 = I(`=cmax') // nw 3
mat define S6 = I(`=cmax') // nw 4 
mat define S7 = I(`=cmax') // no 1
mat define S8 = I(`=cmax') // no 1

***************************
**  5: fill matrix ********
***************************

forvalues k = 1/`=bmax' {
     quietly sum cell if bcell==`k' & bcell != . 
     local i = r(mean)
     forvalues l = 1/`=cmax' {
          quietly sum cell if cell==`l' 
          local j = r(mean)
		  ** run through combinations and put in distance **
		  *mat X[`k',`l']=(x_`i'_`j')
		  mat S1[`k',`l']=(s1_`i'_`j')
		  mat S2[`k',`l']=(s2_`i'_`j')
		  mat S3[`k',`l']=(s3_`i'_`j')
		  mat S4[`k',`l']=(s4_`i'_`j')
		  mat S5[`k',`l']=(s5_`i'_`j')
		  mat S6[`k',`l']=(s6_`i'_`j')
		  mat S7[`k',`l']=(s7_`i'_`j')
		  mat S8[`k',`l']=(s8_`i'_`j')
         }
	 forvalues c = 1/`=cmax' {
	** replace own-cross with zero (identical cell)
	* mat X[`c',`c']=0
	mat S1[`c',`c']=0
	mat S2[`c',`c']=0
	mat S3[`c',`c']=0
	mat S4[`c',`c']=0
	mat S5[`c',`c']=0
	mat S6[`c',`c']=0
	mat S7[`c',`c']=0
	mat S8[`c',`c']=0
        }
    } 

	
***********************************************
****** 6: post results to weight matrix. ******
***********************************************
cd "$store"

	********************
	** Sectors 1-4 *****
	********************

preserve

svmat S1, names(s1_)
svmat S2, names(s2_)
svmat S3, names(s3_)
svmat S4, names(s4_)

save spatial_3_wmat_fan_1234, replace

restore


	********************
	** Sectors 5-8 *****
	********************
preserve

svmat S5, names(s5_)
svmat S6, names(s6_)
svmat S7, names(s7_)
svmat S8, names(s8_)

save spatial_3_wmat_fan_5678, replace

restore

log close 

clear 

exit 